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Vibrational spectra and normal modes of mechanically stable particle packings in three dimensions 
are analyzed over a range of compressions, from near the jamming transition, where the packings lose 
their rigidity, to far above it. At high frequency, the normal modes are localized at all compressions. 
At low frequency, the nature of the modes depends somewhat on compression. At large compressions, 
far from the transition, the lowest-frequency normal modes have some plane-wave character, though 
less than one would expect for a crystalline or isotropic solid. At low compressions near the jamming 
transition, the lowest-frequency modes are neither plane-wave-like nor localized. We characterize 
these differences, highlighting the unusual dispersion behavior that emerges for marginally jammed 
solids. 

PACS numbers: 81.05.Rm, 83.80.Iz, 63.50.+X, 64.70.Pf, 



It is well recognized that the high-frequency vibrations 
in amorphous materials are strikingly different from those 
in crystals. In glasses and other amorphous solids, the 
highest-frequency normal modes are localized in space, 
while in crystals they are extended excitations [l], 0, Q • 
It is also appreciated that even at low frequencies the 
normal modes of disordered systems can be dramatically 
different from the long- wavelength plane waves found in 
ordered materials. The vibrational spectra of amorphous 
solids are characterized by "boson peaks" - extra low- 
frequency modes beyond the long- wavelength plane-wave 
phonons found in crystals. The anomalous modes that 
fall within the boson peak are believed to be responsi- 
ble for the unusual low-temperature properties of glasses, 
such as the plateau in the thermal conductivity [4|. 

Nowhere are the excess low-frequency excitations more 
apparent than in a marginally jammed solid, in which 
a system of particles is compressed to the point where 
they first be gin to touch and form a rigid structure 
0, 0, 0, H SM El G3 Gl- In this system, the den- 
sity of normal modes, instead of vanishing as the fre- 
quency is lowered towards zero, as in a crystal, remains 
constant as shown in Fig. 1(a). This leads to an excess 
in the density of states that diverges at zero frequency. 
This divergence has been interpreted as a boson peak 
Q. Thus, one might expect the marginally jammed solid 
to provide the clearest window into the anomalous low- 
frequency normal modes characteristic of all amorphous 
solids. 

Previously, we have found that the characteristic fre- 
quency and size of the boson peak can be tuned system- 
atically by compressing the marginally jammed solid to 
higher packing fractions Q . In the present paper we an- 
alyze the structure of the normal modes of disordered 
packings in three dimensions (3D) in the marginally 
jammed state and at compressions above this state. The 
system we study here is identical to the one used pre- 
viously to calculate the density of vibrational states 
0, H, 111- We simulate a 3D system of N (1024 < N < 



10000) monodispcrse soft-spheres of mass m and diame- 
ter a in cubic simulation cells employing periodic bound- 
ary conditions. The particles interact via a finite-range, 
purely repulsive, harmonic potential: 
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where r is the center-to-center distance between two par- 
ticles. Length and time are measured in units of a and 
(md 2 /Vq) 1 / 2 respectively. We initially place N parti- 
cles at random in a cubic box of linear dimension L. 
This corresponds to a T = oo configuration. We use 
conjugate-gradient energy minimization [l5[ in order to 
obtain T = configurations. In order to show visu- 
alizations of the normal modes, we have also studied 
N = 10000, bidisperse, 50:50 mixtures of harmonic discs 
with a diameter ratio of 1.4 in 2D. 

The onset of jamming in the limit of large N occurs at a 
packing fraction, </>£° = 0.639±0.003, and is characterized 
by the onset of a nonzero pressure. We determine <j> c for 
each of our finite-system initial configurations by incre- 
mentally compressing (decompressing) until the pressure 
just becomes nonzero (just reaches zero). We then com- 
press the system to obtain zero-temperature compressed 
configurations at controlled values of <j> — <p c . For each of 
these configurations, we compute and diagonalize the dy- 
namical matrix [l6l |. The eigenvalues and eigenmodes of 
this matrix are respectively the squared frequencies, lu 2 , 
of the normal modes of vibration and the corresponding 
polarizations of each particle i in the normal mode 
of frequency u). 

In a previous paper 0, we analyzed the density of 
vibrational states D(u>), of configurations above the jam- 
ming threshold, for systems at <f) > 4>ci an d found three 
characteristic regimes, as labeled in Fig. [TJb). In regime 
A, D(u>) decreases towards zero as u> — > 0. In regime B 
(including B' in Fig. Ufa)), D(u) is approximately con- 
stant, very different than for crystals. Finally, in regime 
C at high frequencies, D(u) decreases with increasing 



frequency. Figure QJa-c) show how the different regimes 
shift with increasing compression. At high A0 = <j) — cj> c , 
regime A extends to fairly high frequencies and regime 
B is small. As the system is decompressed towards the 
marginally jammed solid at A</> = 0, regime A shrinks 
and regime B grows, extending all the way down to zero 
frequency, as indicated by B' in Fig.fjja), while regime C 
remains approximately the same. The growth of regime 
B at the expense of regime A with decreasing compres- 
sion signals the proliferation of anomalous low-frequency 
modes. These are the modes whose structure we wish to 
understand. 
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1: Density of states X>(w) of 3D jammed packings with 
1024 monodisperse spheres at three different compres- 



sions: (a) A0 = 10" 6 , (b) A<f> = 10 , and (c) A<p = 10 . In 
(a) we identify regime B', and in (b) we identify the regimes 
A,B, and C, discussed in the text. 



In order to visualize the nature of the modes, we turn 
to the bidisperse 2D system. Fig. [5] shows typical normal 
modes from regimes A, B', B, and C, defined in Fig.[TJa) 
and (b). In the left panels, the polarization vector for 
each particle is shown, while in the right panel, each par- 
ticle is shaded according to the magnitude of its polar- 
ization vector. The mode from regime A, Fig.[2ftop), ap- 
pears to have some plane-wave-like character, although 
contributions from several different wavevectors are read- 
ily apparent. The high-frequency mode corresponding to 
regime C, Fig. ^bottom) is quite localized. This partic- 
ular mode is representative of the high-frequency modes 
at all values of A<f>; visualizations for different A</> are 
indistinguishable from that shown in Fig. ^bottom) . 

The modes from regimes B' and B, Fig. ^middle), 
are neither plane-wave-like nor localized. The right pan- 
els more clearly reveal the filamentary nature of the ex- 
tended vibrational modes in regimes B' and B. Here we 
point out that these visualizations already suggest some 
subtle differences between regimes B' and B which we 
quantify below (see Fig. [7] and related discussion). 

Fig. [2] suggests that the modes from both regimes A 
and B are extended. To analyze directly how extended 
each mode is, we calculate the inverse participation ratio, 
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FIG. 2: Normal modes of a 2D jammed packing of a 50:50 mix 
of TV = 10000 bidisperse disks, with size ratio 1 : 1.4, interact- 
ing with the potential defined in Eq. [T] These figures corre- 
spond to the regimes identified in Fig. [TJa),(b) characterized 
by(A<£,w). Top row panels: Regime A (lx 10" 1 , 1.76x 10~ 2 ). 
Second row: Regime B' (1 x 10" 6 ,3.35 x 10" 4 ). Third 
row: Regime B (1 x 10" 6 ,0.3). Bottom row: Regime C 
(1 x 10 _1 ,2.30). Left panels: black lines represent the am- 
plitude and direction of the particle vibrations in that mode. 
Right panels: particles are shaded according to the magni- 
tude of their polarization vector. The scale bar indicates that 
darker shading corresponds to a larger ratio of the amplitude 
to the maximum amplitude of particle displacement in that 
mode. Particles with no contact neighbors are not shown. 



mode uj. We show the results for A<f> =10 6 , A<f> = 
1CT 3 , and A</> = lfT 1 in Fig. O We find that on a 
log-log scale, the participation ratio looks quite similar 
for values of A<j> up to at least A(j> sw 10~ 3 . For high 
compressions (A<p = 0.1 and higher), P" 1 at the very 
highest frequencies (regime C) appears to decrease with 
increasing compression, but P~ l is still large in regime 



3 



C at all compressions. Thus, the modes in regime C 
are localized at all compressions. At lower frequencies, 
however, up to lo k, 2, we find P' 1 << 1, indicating that 
the modes are extended over the size of the system. Thus, 
modes in regimes A and B are extended, while those in 
regime C are localized. 
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FIG. 4: Distribution, P(s), of level spacings, s, normalized 
by the the average level spacing, at A(f> = 10 (■) and 10 -6 
(o) in 3D. Inset shows the linear behavior at small s. The line 
is a fit to the Wigner-Dyson function of Eq. [3] 



FIG. 3: Inverse participation ratio, P~ 1 , for 3D jammed pack- 
ings of N = 1024 monodisperse spheres, at three compres- 
sions. 

For these jammed, mechanically stable packings, there 
exist local correlations in the force constants that consti- 
tute the dynamical matrix: The forces around each and 
every particle must be locally balanced. Upon analyzing 
the spacings Alo = u]j + i —uij, between successive normal 
mode frequencies, we find level-repulsion in the distribu- 
tion of level spacings, P(s), where we define s = < ^" > , 
as the level spacing normalized by the average < Alo >. 
The distributions shown in Fig. |4] show little dependence 
on distance to the jamming threshold, and are described 
quite well by the Wigner-Dyson distribution 

P(s) = y e-«> 2 / 4 . (3) 

The distributions shown in Fig. [4] are peaked around the 
average spacing and indicate level repulsion by the lin- 
ear behavior at small spacing. Thus, even though the 
dynamical matrix is sparse, due to the short-range na- 
ture of the interaction potential, the level spacings are 
not completely random, which would lead to a Poisson 
distribution, nor are they completely correlated [l8|. 

Another way to characterize the modes is to look at 
local correlations of the polarizations of neighboring par- 
ticles. We calculate a quantity, cos# e , that is similar to 
the phase quotient parameter often probed in glasses , 

cos 6 e (u) = —— • e ju (4) 

-''pairs ■ - 

where the sum only runs over the number of pairs of 
particles that interact with each other, A pa i rs , and the 



normalized polarization vector of particle i associated 
with the normal mode of frequency lo, = jj^- For 
a mode in which every particle is vibrating in approxi- 
mately the same direction, i.e., strongly correlated mo- 
tion, one would expect cos6* e ss 1. Figure [5] shows cos# e 
as a function of frequency at three values of A(j>. Up to 
compressions of order Acf) = 10~ 3 , the behavior of cos# e 
with lo is insensitive to compression. Over that range 
of compressions, cos 9 e decreases linearly with frequency, 
showing that in the low- frequency modes, particle dis- 
placements are more correlated with their neighbors and 
in high frequency modes, particle displacements are more 
anti-correlated with their neighbors. Note that the fre- 
quency ranges of regimes A and B change appreciably 
with A(j), while cos 9 e (uj) docs not; this suggests that 
the correlations are not noticeably different in the two 
regimes. At high compressions, the curve begins to show 
a kink at approximately lo = 1.75 (the boundary between 
regimes B and C). 

As noted earlier (and by visual inspection of the 2D 
modes of Fig. [2]), the low-frequency modes both near and 
far from the jamming threshold appear to have somewhat 
different character. We have argued that compressed 
systems away from the jamming threshold contain low- 
frequency modes that are more plane- wave-like, whereas 
near the jamming transition, these extended modes are 
very different from plane waves. At intermediate and 
high frequencies the modes appear relatively insensitive 
to packing fraction. In an effort to further differentiate 
between the nature of the low-frequency modes we mea- 
sure the spatial extent of correlated vibrational motions. 
For each mode of frequency lo, we compute the correla- 
tion of polarization vectors between particles i and j, 

C(r ij ) = (e(r i )-e{r j )}. (5) 
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FIG. 5: Correlation between particle displacements in a mode, 
measured using cos 9 e denned by Eq. [4] at three compres- 
sions for 3D jammed packings with TV = 1024 monodisperse 
spheres. 



FIG. 7: Dependence of the characteristic length scale £, de- 
fined in text, on the normal mode frequency u>, at three values 
of A(f> = 1CT 6 (O), 10" 3 (□), and KT 1 (0), for N = 10000 
monodisperse spheres in 3D. 



In Fig. [S] we show C(r) for the lowest lying frequency 
modes at three different compressions, A(f> = 10 -6 , A<fi = 
10 -2 , and Acf> = 10 _1 , for our 3D jammed packings. 



(b) 




FIG. 6: Spatial correlation function C(r) of particle vibrations 
at three values of A<t> = 10~ 6 (solid line), 10 -2 (dashed), and 
10 _1 (dotted), for TV = 10000 monodisperse spheres in 3D. 
(a) Lowest-lying frequencies, and (b) intermediate frequen- 



For the low-frequency modes we find somewhat 
stronger correlations at higher compressions. For the 
middle-to-high frequency range of the vibrational spec- 
trum, beyond regime A, the normal modes become in- 
creasingly similar at different compressions. In regime 
C, the modes are indistinguishable at different compres- 
sions. 

All of the correlation functions in Fig. [5] decay non- 
monotonically to zero and cross zero at some finite r. 
We define the value of r at which C(r) first crosses zero 
to be t(us). In a pure plane wave, particle vibrations are 
correlated and C{r) will cross zero at the scale of half 
the wavelength, so for ordinary sound modes at low fre- 
quency, one would expect l(u>) oc l/u>. In Fig. [Jj we 
plot I as a function of to for different A<\> at low-to- 
intermediate frequencies. For the system closest to the 
jamming threshold, t is approximately independent of uj 



at very low frequencies. Beyond this constant region, 
I decreases with increasing frequency, corresponding to 
moving along the plateau in the density of states from 
regime B' to B [20(1. As the system moves further from 
the jamming threshold, i.e., as A<j) increases, the region of 
constant I shrinks. Over the range of frequencies where I 
decreases, the characteristic length is greater the further 
the system is from the jamming threshold. This is to be 
expected as the modes contain more plane-wave charac- 
ter at larger compressions. At slightly higher frequen- 
cies still, the curves begin to overlap, so that the modes 
do indeed become practically indistinguishable. The fre- 
quency at which this occurs coincides with the point in 
the density of states where the plateau regions start to 
merge (see Fig. [T] of Ref. 0). These data suggest that 
the distinction between modes from regimes B' and B is 
related to how extended the modes are. This is already 
evident from the visualizations presented in Fig. [2] 

Another way to quantify the difference between the ex- 
tended modes in regimes A and B is through the Fourier 
transforms of the eigenmodes at different frequencies 
uj throughout the spectrum. Specifically, we take the 
Fourier transform of the appropriate component, either 
longitudinal or transverse, of the particle polarization 
vector ej(uj), of each particle j [1, @|: 



(6) 



h{k,u) = (j? XVk-e^exp(ik-rj) 
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k A ej u cxp(ik • rj ) 



In a perfect crystal, these functions would be delta- 
functions at the wavevectors k in each Brillouin zone 
characterizing the longitudinal or transverse vibrational 
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modes at frequency u>. In Fig. [8j we show fi,(Jk, u>) and 
/^(fc, lo) curves for two disordered configurations in three 
frequency bands: (i) at the lowest frequency, (ii) in the 
middle of the band, and (iii) at the high-frequency end 
of the spectrum. Each curve is averaged over a nar- 
row bin of frequencies. For comparison, Figs. [8jja),(b) 
are from a system very close to the jamming thresh- 
old, at Acj) = 10 -6 , and Fig. [^c),(d) are for a sys- 
tem that is highly compressed and far from the jamming 
threshold, at Ac/> = 10 _1 . At A</> = 1CT 6 , the low and 
mid-frequency curves correspond to vibrational states in 
regime B' defined in Fig.QJa), where the density of states 
is relatively flat, while the high frequency curve corre- 
sponds to regime C. At high compression, A<fr = 10 _1 , 
the low frequency curve corresponds to a state in regime 
A, the mid-frequency curve corresponds to regime B and 
the high frequency curve corresponds to regime C. The 
longitudinal functions in general show much more pro- 
nounced structure than do their transverse counterparts. 
The only exception to this occurs at very low frequencies 
and small wavevectors. In this region the transverse func- 
tion has a very tall first peak and the longitudinal func- 
tion shows only a hint of structure. The low-wavevector 
part of the peak in fT(k,LU\ ow ) is not resolved because 
it occurs at k < 2ir/L where L is the size of the system. 
That is, the peak is cut off because of the finite size of the 
system. The first peak in /i(fc,o;iow) is also absent for 
the same reason. In order to see this structure, one would 
have to either use a larger box size at the same value of 
frequency or else look at uj) at a larger value of w. 



3 




0.05 



0.1 



3 



0.05 







— Low-co 






--- Middle-<o 






— High— co 










H.J ~" w ( C ) 



3 



0.1 



0.05 




5 10 15 20 

k 



5 10 15 20 

k 



FIG. 8: Fourier transforms of the eigenmodes for the low 
(solid line), middle (dashed line), and high frequency (dot- 
ted line) regions of the vibrational spectrum, at two extreme 
compressions in 3D. Top panels: A(f> = 1 x 1CP 6 , for the, 
(a) longitudinal, and (b) the transverse components. Bottom 
panels: A<f> = 1 x 10 , for the, (c) longitudinal, and (d) 
the transverse components. Insets to (b) and (d) show the 
dominance of the low-fc peak in the transverse functions. 

There are multiple oscillations visible in the longitu- 
dinal response, /i(fc,w). This structure can be thought 



of as the equivalent of the repeated structure seen in the 
higher Brillouin zones of a crystal 0]. It reflects the large, 
sharp first peak in the pair-correlation function, g(r) [211 ] . 
which leads to strong oscillations in the structure factor 
= (ISjCxP^k ' r )| 2 ) 0- Similar oscillations show 
up but with a much smaller amplitude in fT(k,u>). 

Overall, the results look fairly similar for the two com- 
pressions. For the longitudinal response, the peaks are 
somewhat smaller and wider at low compression. How- 
ever, the most obvious difference is not in the peaks but 
in the minima between them which become more shallow 
at small A</>. This is particularly apparent at the lower 
frequencies in the longitudinal response. This means that 
more wavevectors are making significant contributions to 
the low-frequency longitudinal modes at low compression 
(which arc in regime B) than to low-frequency longitu- 
dinal modes at high compression (which are in regime 
A), making the mode very different from any plane wave 
with a single wavevector. A similar trend is apparent 
in the transverse response. Contributions of wavevec- 
tors different from the peak value are relatively larger at 
low compressions; that is, more wavevectors contribute to 
eigenmodes in regime B than in regime A. The intermedi- 
ate wavevector oscillations, clearly visible at A<fi = 10 _1 , 
are much less pronounced near the onset of jamming. At 
low frequency, the transverse response is practically flat 
at wavevectors k > 5 - all of these high wavevectors con- 
tribute nearly equally in regime B. 

The velocity of longitudinal or transverse sound can 
be estimated from the frequency-dependence of the po- 
sition of the first peak in _/x(fc,w) or fr{k, oS). Unlike in 
crystals, the dispersion curve is not well-defined because, 
as we have seen in Fig. [FJ w) and fj-{k,uS) have 

significant amplitude over a large range of wavevector at 
all frequencies. As a result, it is not sufficient to simply 
plot the positions of the peaks of the dynamic structure 
factors as one might do for a crystal. 

To emphasize the idea of mode mixing, the dispersion 
data are best visualized on a gray scale plot as shown 
in Fig. [9] (The features outlined here have been seen 
in a ran ge o f glassy systems, including soft-sphere [12] , 
covalent [23||, and metallic [24j glasses.) Because the am- 
plitude of the first peak in ff(k,u)) is so strong, varia- 
tions at wavevectors greater than the first peak are lost 
in Fig. O Therefore, to observe the underlying structure 
in the dispersion data at larger wavevectors, in Fig. 1101 
we show the transverse dispersion curves over a limited 
range in amplitude. These nicely contrast the unusual 
underlying structure to the dispersion relations at the 
two compressions. 

It is clear from looking at these plots that the disper- 
sion relations are very broad indeed, especially at low 
compression. Moreover, there is very little difference be- 
tween the two compressions, although the contrast de- 
creases for both transverse and longitudinal modes with 
decreasing compression. The limit of what we can achieve 
is A(j) = 10~ 6 . It is not clear to us if one were able to go 
even closer to A<j) — whether the variations in /r(fc,u;) 
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FIG. 9: Dispersion relation at two different compressions in 
3D. Left panels: A</> = 1 x 1(T 6 . Right panels: A(f> = lx 1CT 1 . 
Top and bottom panels correspond to longitudinal and trans- 
verse modes respectively. Darker shading corresponds to 
larger amplitude in /l,t- The maximum of the transverse 
amplitudes are typically a factor of 5 larger than the longi- 
tudinal data. The bars to the right indicate the amplitude 
scale. 




FIG. 10: Transverse dispersion data in 3D for A(j> = 1 x 10 
(left) and 1 x 10" 1 (right), over a limited range in amplitude 
/t (see scale bar). 

would disappear entirely. 

It is also unclear how to define a proper velocity of 
sound not only because the peaks in the dispersion re- 
lations are so broad, but also because their amplitudes 
decay rapidly with increasing (k,u>). We illustrate this 
latter point in Fig. QT] where we show the approximately 
exponential decay of the maximum peak height of the 
transverse Fourier modes with increasing frequency. We 
have been unable to determine precisely whether the de- 
cay constant depends on compression; although there ap- 
pears to be a small difference between the two compres- 
sions shown in Fig. [TTJ that difference is small. 

In summary, we have studied the characteristics of vi- 
brational modes in different frequency regimes. From 



the density of states, we see that there are three regimes: 
regime A, where the density of states drops towards zero 
with vanishing frequency; regime B, where the density of 
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FIG. 11: Decay of the peak maximum of the transverse mode 
transforms within the first "Brillouin zone", for A(f> = 10 _1 
(solid line) and 10 -6 (dashed line) in 3D. 



states is approximately flat; and regime C, where the den- 
sity of states is decreasing towards zero with increasing 
frequency. As the system is decompressed towards the 
marginally jammed state, regime B (or B') increases at 
the expense of regime A while regime C is relatively un- 
affected. Modes from regime C arc localized while modes 
from regime A and B arc extended. Those in Regime A 
are somewhat more plane- wave- like, in that contributions 
from wavevectors away from the peak of the dynamic 
structure factor are fairly small (the peak is relatively 
narrow) while those in Regime B have broader peaks 
in the dynamic structure factors, with significant con- 
tributions from a wider range of wavevectors. Generally, 
however, we do not observe strong differences between 
modes from regime B and regime A. Thus, the change 
in the nature of the modes with decreasing compression 
is much less dramatic than the change in the density of 
vibrational states. 
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